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We calculate the spectrum of ultra high energy cosmic rays produced by 
the decay of a superheavy dark matter population clustered in the galactic 
halo. To perform this calculation we start with fragmentation functions mea- 
sured at LEP and evolve them to the cosmic ray energy scale using the QCD 
DGLAP equations. We consider Standard Model evolution and supersymmet- 
ric evolution. We also take into account many-body final states in the decay 
of the dark matter particles. 
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1 Motivation 



Ultra High Energy Cosmic Rays (UHECR) are microscopic particles -protons, photons or 
perhaps more exotic objects- with a macroscopic energy, about 50 Joules per particle in 
the extreme of the present observed spectrum. They strike the upper layers of the Earth 
atmosphere at a rate of about one event per century and per kilometer square. Over one 
hundred of them have been detected so far by kilometer scale detectors and many more 
are expected to be seen by forthcoming observatories. 

Explaining how these particles get this huge energy is a real challenge for our present 
understanding of the cosmos. One UHECR primary gets somehow ten million times the 
energy that a proton will gain at the future Large Hadron Collider (LHC) at CERN. 
Nature has always been at ease beating the achievements of mankind. One can envisage 
two broad classes of mechanisms by which nature can produce UHECRs, the so-called 
bottom up class models and the top-down class models. 

In bottom-up models charged particles are accelerated by magnetic fields in large 
astrophysical sites. In top-down models particles are not accelerated but are created at 
birth with the huge energy typical of UHECRs. 

Here we will concentrate on a particular top-down model |], |2[ |3[] . We will assume the 
UHECRs are produced by the decay of a population of superheavy dark matter particles 
with a lifetime longer than the age of the Universe r > 10 10 yr and with mass ~ 100 
Joule/c 2 . Theoretical motivation for these particles can be found in |4]]. 

A superheavy dark matter population created at some stage of the early universe || |(| 
will gravitationally cluster in the galactic halo. Since the length scale of the halo is around 
100 kpc, UHECRs produced in the galactic halo will not have time to interact with the 
CMB before they reach the Earth. The absence of GZK cut-off is a genuine prediction of 
models where UHECRs are produced by the decay of a superheavy dark matter particle 
clustered in the galactic halo. 

Since the halo shape is close to spherical one expects a quasi isotropic distribution of 
events from a halo superheavy dark matter population, which is compatible with experi- 
ments. At present it is not possible to make any strong claim about the observed angular 
distribution of UHECR events because of low statistics. Hopefully, future observatories 
like Pierre Auger will gather a large enough sample of events to settle down this issue |7J. 

Finally there is the question of the UHECR composition. The cosmic ray observatories 
cannot measure the composition of the primary flux (whether they are photons, protons 
or heavy nuclei) in an event by event basis. Composition can only be determined in a 
statistical way. Present analyses tend to favour protons as the main component in the 
primary flux ||. In top-down models the main component in the production site, in 
our case the the galactic halo, is the photon component (actually, neutrinos dominate 
over photons in a large range of the spectrum but their probability to be detected is too 
small with present or past detectors). However, on their way to the Earth, photons will 
interact with the low frequency radio background in the Galaxy and their total flux may 
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be substantially diminished so that on the Earth the total photon flux may be comparable 
or smaller than the baryon flux. Whether this is indeed possible, subject to the EGRET 
bound on the low energy 7-rays which will be created by the electromagnetic cascading, 
will be discussed elsewhere. 

Summing up, there are three main tests to falsify or support the hypothesis of UHECRs 
produced by the decay of dark matter particles clustered in the galactic halo. The first one 
is the energy distribution of events or spectrum. The second one is the expected angular 
distribution of events in the sky. The third one is the composition of the primary flux. 
Here we will focus on the first test and, partially, on the third one. We will briefly show how 
quantum chromodynamics (QCD) can be used to calculate the spectrum and composition 
(without photon galactic processing) of the expected UHECR flux. For further details 



see pL M. QCD was also used to calculate the spectra of UHECRs in [11, 12, 13 



2 DGLAP Evolution 

Let us assume that a dark matter particle X population, superheavy (M x > 3 x 10 12 GeV) 
and metastable, with a lifetime larger than the age of the universe, is clustered in the 
galactic halo. Barring unnatural relations between the coupling constant of X and the 
other fields the total particle multiplicity produced by its decay will be dominated by 
partonic decays. The flux for the primary h (baryon, photon or neutrino) is proportional 
to the inclusive decay width 



1 dT(X^h + ...) _ fidz 1 dT a (y^\M 2 x ) 



J h (E)oc-^ y - ' = V / — 

T x dx V Jx z r a dy 



y=x/z 



Particle h carries a fraction x of the maximum available momentum Mx/2, and a fraction 
z of the parton a momentum. The first term in the integrand is the decay width of X 
into parton a, dT a /dy, which is calculable in perturbation theory; in lowest order and for 
2-body decay it is proportional to 5(1 — y)- The second factor, the non-perturbative D^, 
is the fragmentation function (FF) for particles of type h from partons of type a. It gives 
the expected mean number of particles h coming from parton a. Fragmentation functions 
cannot be calculated from firsts principles but their dependence in the energy scale fi is 



governed by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [|14|. [15 



dD ^P = E «.0« a )) ® ^ ( 2 ) 

where a s (/i 2 ) is the strong coupling constant and Pb a (x,a s ) is the splitting function for 
the parton branching a — > b. Here the convolution of two functions A(x) and B(x) is 
defined as 

A(x) ® B(x) = / —A{z)B{-). (3) 
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The splitting functions can be expanded perturbatively: 



P ba (x, a s ) = P ba {x) + 0(a s 



(4) 



We limit our study to leading order in a s and therefore ignore 0(a s ) corrections to the 
splitting functions. It is also convenient to define the following dimensionless evolution 
parameter 
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2ixb a s (fi 2 ) 

b being the coefficient in the leading order /3-function governing the running of the strong 
coupling: /3(a s ) = —ba 2 s . We take D\ to represent the sum of particle h and, if different, 
its antiparticle h 

In the Standard Model there are two partons species: gluons g and quarks q (for the 
sake of brevity we shall only consider the singlet quark, i.e., the sum of all quark and 
antiquark flavours). The DGLAP equations can be written as 
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When supersymmetry (SUSY) is included one has in addition the gluinos and the 
squark singlet. The DGLAP equations are now 
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We have developed a C++ code to solve the DGLAP equations for fragmentation func- 
tions. We have extended the algorithm introduced in Ref. [0 to be able to solve the 
SUSY equations as well as the Standard Model equations [[K] . 

We begin with FFs at the energy scale fiQ = Mz and evolve them to the final energy 
scale fi = Mx using Eq. (|6]) and also Eq (0) if SUSY is included. For the initial fragmenta- 



tion function of baryons, D%(x, Mf ) + D^(x, Mf ), we adopt the fit performed in Ref. [|Tl 
to LEP hadronic data JT7|]. For photons and neutrinos we generate initial data at the Z 





peak using the QCD Monte Carlo event generator HERWIG 



data shows that although HERWIG overproduces baryons at high x |TT| , |1 



Comparison with LEP 
its photon 



and meson output at the Z peak matches the experimental spectra remarkably well. Since 
neutrinos mainly come from charged pion and kaon decays, one can thus be confident in 
taking the HERWIG generated FF for neutrinos as the initial condition for the evolution. 
There is also a sizable contribution from heavy flavour decays to the neutrino spectrum 
at high x which is explicitly taken into account by HERWIG. 

A (s)parton is not included in the evolution as long as the energy scale is lower than 
its mass; when the threshold for its production is crossed, it is added to the evolution 
equations with an initially vanishing FF and it is assumed to be a relativistic particle. 
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In the SM case we evolve the q and g initial fragmentation functions from Mz to M t , 
the top quark mass, with the number of flavours set to np = 5, and then evolve from M t 
to Mx with n-p = 6. 

In the SUSY case we evolve the q and g initial fragmentation functions from Mz 
to the supersymmetry breaking scale Msusy > M t using the SM equations to obtain 
D i( x > m susy)> with i = q,g. Then we take D*(x, M| USY ), i = q,g, and D}(x, M| USY ) = 0, 
j = s, A, and evolve them from MgusY to M x using the SUSY equations. All spartons 
are taken to be degenerate with a common mass Msusy- 



3 Ultra High Energy Cosmic Ray Spectrum 

We can now translate the calculated fragmentation functions into the expected cosmic 
ray spectrum in order to confront the observational data. In the previous Section (HD we 
have briefly shown how to calculate quark singlet and gluon functions for the SM, and 
quark singlet, gluon, squark singlet and gluino functions for SUSY. In the absence of a 
specific model for the different branching ratios we weight all the (s)parton contributions 
evenly. For 2-body decay of X, x = 2E/Mx the flux of particle h is given by ||: 

E 3 J halo (£) = Bx 3 D h (x, M 2 X ). (8) 

We have multiplied the flux by E 3 , as is usual, to emphasise the structure in the spectrum 
near the GZK energy. The normalisation factor B is common for the galactic halo flux of 
baryons, neutrinos and photons, and determines the quantity ux/tx (see Eq. |l|). 

Let us now compare the calculated cosmic ray flux to the published data from Fly's 
Eye |2(|, AGASA 0, Haverah Park [|2| and Yakutsk |23|. In Ref. @ these data have 



been carefully assessed for mutual consistency and appropriate adjustments made to the 
energy calibration. Its authors recommend adoption of the following standard differential 
energy spectrum below the GZK energy, in the range 4 x 10 17 eV < E < 6.3 x 10 18 eV: 

/ E \ -3.20±0.05 

J(E) = (9.23 ± 0.65) x 1(T 33 m~ W^V" 1 - - ^ - j , (9) 

with the spectrum flattening at higher energies as J(E) oc E~ 2 ' 75±0 ' 2 up to the GZK 
energy, and extending further to at least 3 x 10 20 eV |24j]. Thus the UHECR spectrum 
can naturally be interpreted |2D| as the superposition of the 'low energy' component (0), 
and the new 'flat' component that extends into the post-GZK region. The former is 
presumably galactic in origin (consistent with the detection of anisotropy at ~ 10 18 eV 
|25| , |26|j), while the latter is interpreted as produced by the decay of a superheavy 
particle population in the galactic halo. Taking baryons to be the dominant primary 
UHECRs as indicated by experiment, the total flux is 

E 3 J(E) = A + Bx 3 D b ^(x, M|), (10) 
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where the values of k and m can be read off Eq. (||). Note that since £) bar y° n and D 1 have 
a similar shape, taking photons to be the primaries would just alter the normalisation B. 
In Fig. [I] we plot the best SM evolution fit to the cosmic ray data while in Fig. |] we plot 
the best SUSY evolution fit. 
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Figure 1: The best SM evolution fit to the cosmic ray data with a decaying particle mass 
of 10 12 GeV. The dotted line indicates the extrapolation of the power-law component from 
lower energies, while the dashed line shows the decay spectrum; the solid line is their sum. 



The assumption of 2-body decay may be rather naive for a superheavy particle like 
a crypton Q which is expected to decay through very high-order non-renormalisable 
operators. Many-body decay distributes the total energy Mx among several particles and 
thus flattens the spectrum. We assume that many-body effects are purely kinematical and 
hence can be encapsulated in the phase space of the decay. Let p n (z) be the probability 
density that one parton carries off a fraction z of the total available energy per parton 
Mx/2. For n > 3 we get 

p n (z) oc z{l - z) n -\ (11) 
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Figure 2: The best SUSY evolution fit to the cosmic ray data with a decaying particle 
mass of 5 x 10 12 GeV and M$usy — 400 GeV. The dotted line indicates the extrapolation 
of the power-law component from lower energies, while the dashed line shows the decay 
spectrum; the solid line is their sum. 



To leading order in a s the particle flux is then given by 

dz 



E 3 J h ^°(E)=Bx 3 [ —p n (-)D h (z,M 2 x ). 

Jx z V z J 



(12) 



In particular if D h (x, M|) oc (1 — x) a ( M x) as x — > 1, the differential particle flux decreases 
as J halo (x) oc (1 - x )«(a4)+™-2. 

For many-body X decays the total flux is 



■>^z,M 2 x ), 



(13) 



where n is the number of partons into which X decays. In Fig. |] we plot evolved SUSY 
spectra for different values of n. 
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Figure 3: Cosmic ray data compared with SUSY evolved spectra for a decaying particle 
mass of 10 13 GeV with M S usy = 400 GeV, and many-body decays to n partons: n = 2 
(solid line), n = 8 (dotted line) and n — 16 (dashed line). 



4 Conclusions 



We have calculated the UHECR spectra of baryons, photons and neutrinos expected 
from the decay of superheavy dark matter particles of mass Mx ~ 10 12 GeV which are 
clustered in the galactic halo. The spectrum for every primary particle is given by FFs 
at the scale Mx- We have calculated these FFs using QCD DGLAP equations. We 
have considered Standard Model equations and SUSY equations. We have also taken 
into account the possibility of many-body decay in the decay of X. The shape of the 
fragmentation spectrum (of either baryons or photons) fits rather well the new component 
of ultra-high energy cosmic rays extending beyond the GZK energy. 
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